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Foreward 


I am delighted to introduce the first book on Multimedia Data Mining. When 
I came to know about this book project undertaken by two of the most active 
young researchers in the held, I was pleased that this book is coming in early 
stage of a held that will need it more than most helds do. In most emerging 
research helds, a book can play a significant role in bringing some maturity to 
the held. Research helds advance through research papers. In research papers, 
however, only a limited perspective could be provided about the held, its 
application potential, and the techniques required and already developed in 
the held. A book gives such a chance. I liked the idea that there will be a book 
that will try to unify the held by bringing in disparate topics already available 
in several papers that are not easy to find and understand. I was supportive 
of this book project even before I had seen any material on it. The project 
was a brilliant and a bold idea by two active researchers. Now that I have it 
on my screen, it appears to be even a better idea. 

Multimedia started gaining recognition in 1990s as a held. Processing, 
storage, communication, and capture and display technologies had advanced 
enough that researchers and technologists started building approaches to com- 
bine information in multiple types of signals such as audio, images, video, and 
text. Multimedia computing and communication techniques recognize corre- 
lated information in multiple sources as well as insufficiency of information in 
any individual source. By properly selecting sources to provide complemen- 
tary information, such systems aspire, much like human perception system, 
to create a holistic picture of a situation using only partial information from 
separate sources. 

Data mining is a direct outgrowth of progress in data storage and process- 
ing speeds. When it became possible to store large volume of data and run 
different statistical computations to explore all possible and even unlikely cor- 
relations among data, the field of data mining was born. Data mining allowed 
people to hypothesize relationships among data entities and explore support 
for those. This field has been put to applications in many diverse domains and 
keeps getting more applications. In fact many new fields are direct outgrowth 
of data mining and it is likely to become a powerful computational tool. 




Preface 


Approximately 17 million people in the USA (6% of the population) and 140 
million people worldwide (this number is expected to rise to almost 300 mil- 
lion by the year 2025) suffer from diabetes mellitus. Currently, there a few 
dozens of commercialised devices for detecting blood glucose levels [1], How- 
ever, most of them are invasive. The development of a noninvasive method 
would considerably improve the quality of life for diabetic patients, facilitate 
their compliance for glucose monitoring, and reduce complications and mor- 
tality associated with this disease. Noninvasive and continuous monitoring of 
glucose concentration in blood and tissues is one of the most challenging and 
exciting applications of optics in medicine. The major difficulty in develop- 
ment and clinical application of optical noninvasive blood glucose sensors is 
associated with very low signal produced by glucose molecules. This results in 
low sensitivity and specificity of glucose monitoring by optical methods and 
needs a lot of efforts to overcome this difficulty. 

A wide range of optical technologies have been designed in attempts to 
develop robust noninvasive methods for glucose sensing. The methods include 
infrared absorption, near-infrared scattering, Raman, fluorescent, and thermal 
gradient spectroscopies, as well as polarimetric, polarization heterodyning, 
photonic crystal, optoacoustic, optothermal, and optical coherence tomogra- 
phy (OCT) techniques [1-31]. 

For example, the polarimetric quantification of glucose is based on the 
phenomenon of optical rotatory dispersion, whereby a chiral molecule in an 
aqueous solution rotates the plane of linearly polarized light passing through 
the solution. The angle of rotation depends linearly on the concentration of the 
chiral species, the pathlength through the sample, and the molecule specific 
rotation. However, polarization sensitive optical technique makes it difficult 
to measure in vivo glucose concentration in blood through the skin because of 
the strong light scattering which causes light depolarization. For this reason, 
the anterior chamber of the eye has been suggested as a sight well suited for 
polarimetric measurements, since scattering in the eye is generally very low 
compared to that in other tissues, and a high correlation exists between the 
glucose in the blood and in the aqueous humor. The high accuracy of anterior 
eye chamber measurements is also due to the low concentration of optically 
active aqueous proteins within the aqueous humor. 

On the other hand, the concept of noninvasive blood glucose sensing using 
the scattering properties of blood and tissues as an alternative to spectral 
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absorption and polarization methods for monitoring of physiological glucose 
concentrations in diabetic patients has been under intensive discussion for 
the last decade. Many of the considered effects, such as changing of the size, 
refractive index, packing, and aggregation of RBC under glucose variation, are 
important for glucose monitoring in diabetic patients. Indeed, at physiological 
concentrations of glucose, ranging from 40 to 400 mg/dl, the role of some of the 
effects may be modified, and some other effects, such as glucose penetration 
inside the RBC and the followed hemoglobin glycation, may be important 
[30-32], 

Noninvasive determination of glucose was attempted using light scattering 
of skin tissue components measured by a spatially-resolved diffuse reflectance 
or NIR frequency-domain reflectance techniques. Both approaches are based 
on change in glucose concentration, which affects the refractive index mis- 
match between the interstitial fluid and tissue fibers, and hence reduces scat- 
tering coefficient. A glucose clamp experiment showed that reduced scattering 
coefficient measured in the visible range qualitatively tracked changes in blood 
glucose concentration for the volunteer with diabetes studied. 
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1.1 Introduction 

The Integrated Health Management (IHM) for the future aerospace sys- 
tems requires to interface models of multiple subsystems in an efficient and 
accurate information environment at the earlier stages of system design. The 
complexity of modern aeronautic and aircraft systems (including e.g. the 
power distribution, flight control, solid and liquid motors) dictates employ- 
ment of hybrid models and high-level reasoners for analysing mixed contin- 
uous and discrete information flow involving multiple modes of operation in 
uncertain environments, unknown state variables, heterogeneous software and 
hardware components. 

To provide the information link between key design/performance parame- 
ters and high-level reasoners we rely on development of multi-physics perfor- 
mance models, distributed sensors networks, and fault diagnostic and prog- 
nostic (FD&P) technologies [38] in close collaboration with system designers. 
The main challenges of our research are related to the in-flight assessment of 
the structural stability, engine performance, and trajectory control. The main 
goal is to develop an intelligent IHM that not only enhances components and 
system reliability, but also provides a post-flight feedback helping to optimize 
design of the next generation of aerospace systems. Our efforts are concen- 
trated on several directions of the research. One of the key components of 
our strategy is an innovative approach to the diagnostics/prognostics based 
on the real time dynamical inference (DI) technologies extended to encom- 
pass hybrid systems with hidden state trajectories. The major investments 
are into the multiphysics performance modelling that provides an access of 
the FD&P technologies to the main performance parameters of e.g. solid and 
liquid rocket motors and composite materials of the nozzle and case. 

Some of the recent results of our research are discussed in this chapter. We 
begin by introducing the problem of dynamical inference of stochastic non- 
linear models and reviewing earlier results. Next, we present our analytical 
approach to the solution of this problem based on the path integral formu- 
lation. The resulting algorithm does not require an extensive global search 
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for the model parameters, provides optimal compensation for the effects of 
dynamical noise, and is robust for a broad range of dynamical models. In the 
following Section the strengths of the algorithm are illustrated illustrated by 
inferring the parameters of the stochastic Lorenz system and comparing the 
results with those of earlier research. Next, we discuss a number of recent 
results in application to the development of the IHM for aerospace system. 

Firstly, we apply dynamical inference approach to a solution of classical 
three tank problems with mixed unknown continuous and binary parameters. 
The problem is considered in the context of ground support system for filling 
fuel tanks of liquid rocket motors. It is shown that the DI algorithm is well 
suited for successful solution of a hybrid version of this benchmark problem 
even in the presence of additional periodic and stochastic perturbation of 
unknown strength. 

Secondly, we illustrate our approach by its application to an analysis of the 
nozzle fault in a solid rocket motor (SRM) . The internal ballistics of the SRM 
is modelled as a set of one-dimensional partial differential equations coupled to 
the dynamics of the propellant regression. In this example we are specifically 
focussed on the inference of discrete and continuous parameters of the nozzle 
blocking fault and on the possibility of an application of the DI algorithm to 
reducing the probability of “misses” of an on-board FD&P for SRM. 

In the next section re-contact problem caused by first stage/upper stage 
separation failure is discussed. The reaction forces imposed on the nozzle of the 
upper stage during the re-contact and their connection to the nozzle damage 
and to the thrust vector control (TVC) signal are obtained. It is shown that 
transient impact induced torque can be modelled as a response of an effective 
damped oscillator. A possible application of the DI algorithm to the inference 
of damage parameters and predicting fault dynamics ahead of time using the 
actuator signal is discussed. 

Finally, we formulate Baesian inferential framework for development of 
the IHM system for in-flight structural health monitoring (SHM) of composite 
materials. We consider the signal generated by piezoelectric actuator mounted 
on composite structure generating elastic waves in it. The signal received by 
the sensor is than compared with the baseline signal. The possibility of damage 
inference is discussed in the context of development of the SHM. 


1.2 Dynamical inference of stochastic nonlinear models 

Complex phenomena in nature and technology can often be modeled suc- 
cessfully by stochastic nonlinear dynamical systems, thereby facilitating the 
diagnosis of faults, the prognosis of future conditions, and control. Exam- 
ples range from models from molecular motors [47] to epidemiology [18] and 
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from solid rocket motors [40] to coupled matter radiation systems in astro- 
physics [14]. 

The problem of inferring the parameters of such models from time-series 
data has therefore attracted much attention over the last decade. Although 
no general method exists for inferring the parameters of stochastic nonlinear 
dynamical models from measurements, various schemes have been proposed 
[33, 20, 22, 35, 19, 42, 48, 45]. An important numerical technique, suggested 
in [43, 20], is based on estimating drift and diffusion coefficients at a number 
of points in the phase space of the dynamical system. A particle filter ap- 
proach [23] and the Markov Chain Monte-Carlo (MCMC) approach [13] were 
applied successfully to reconstruct the model parameters alone. 

Arguably the most general approach to the solution of this problem can 
be obtained within Bayesian inferential framework [22, 35, 11]. The Bayesian 
method was used for parameter estimation in maps in the presence of dy- 
namical [34] and weak measurement [35] noise and was applied to inference of 
continuous systems in [19]. 

A common drawback of these earlier works is their reliance on numerical 
methods for the optimization of cost functions. This disadvantage becomes in- 
creasingly more pronounced when inferring on-board complex aerospace sys- 
tems. Another deficiency is that most of the earlier works deal with discrete 
maps, and the corresponding results are therefore not immediately applicable 
to continuous systems. This problem was addressed in [19], however, an ad 
hoc likelihood function was used there. 

We now formulate the problem and present its solution using path integral 
approach within Bayesian inferential framework. 

The time-series data of the control variables is usually observed at sequen- 
tial time instants and are related to the unknown “tru” system L-dimensional 
states via the measurement model. In these settings the following formula- 
tion of the problem of dynamical inference is commonly adopted. Given M- 
dimensional time-series data y = {y„ = y(t n )} (t n = nh), how can one 
infer the time variation of the unknown model parameters and the unknown 
dynamical trajectory A4 = |c(t), b(t), D, M, {x„}|? It is assumed that the 
underlying dynamics can be described by a set of L-dimensional (L > M ) 
stochastic differential equations of form 

x(t) = f(x\c) + VD£(t), (1.1) 

and the observations y are related to the actual unknown dynamical variables 
X = {x„ = x(f„)} via the following measurement equation 

y(t) = g(x|b) + VM?7(t). (1.2) 

Here X is an M x L measurement matrix, £(f) and r/(f) are L- and M- 
dimensional Gaussian white noises, and D and M are L x L and M x M 
dimensional dynamical and measurement diffusion matrices respectively. 
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The solution of this problem is given by the so-called posterior probability 
density function (PDF) p p0 st (AT |3^) of the unknown parameters A4 condi- 
tioned on observations. In Bayesian model inference, the posterior PDF is 
related to the prior p pr {M) PDF via Bayes’ theorem [15]: 


P pos t{M\y) 


e(y\M) Ppiiol fM) 
f i{y\M) p prio r(M)dM 


(1.3) 


Here, the likelihood function £(y\M) is the probability density to observe 
{y„(t)} given choice A4 of the dynamical model. Meanwhile, the prior acts 
as a regularizer, concentrating the parameter search to those regions of the 
model space favored by our expertise and any available auxiliary information. 
This initial assignment of probabilities must be “coherent” [24], i.e., consistent 
with the physics of the problem. In practice, (1.3) can be applied iteratively 
using a sequence of data blocks y, y ' , . . .; the posterior computed from block 
y serves as the prior for the next block y ' , etc. For a sufficiently large num- 
ber of observations, p ps {M\y ,y' , . . .) becomes sharply peaked around a most 
probable model M* . 

If the sampling is dense enough the problem can be conveniently solved 
using Euler mid-point discretization of Eqs. (1.1), (1.2) in the form 


x„+i = x„ + hf(x*Jc) + Vhi>t n , 
y n = g(x n |b) + Mt 7 „, 


where x* = (x n+ i + x„)/2. It was shown earlier (see e.g. [?, 46]) that for 
independent sources of white Gaussian noise in (1.4) the probability to observe 
y n at each time step can be factorized and written in the form 

0 (v . , l-v ,•) — f dx "~ 1 p-My~+i-g(x„+i|b)] T M - 1 [y„ + i-g(x„ + i|b)l 

^y n+ 1 |x n ,cj J e 

x 1 e -| [in-f(x* |c)] T D _1 [i n -f(x*|c)]-|V.(f(x„)|c) n 

V / (277/ l )i|D| 

Summation over all the discretization points n = 0...N — 1 yields the following 
result for the minus log-likelihood function S = Sd yn + S meas = — \nl{y\M) 


S = 


+ 



h N ~ 1 

2 5Z { V -( f ( x «)l c ) 

n = 0 


[x„ — f(x* |c)] T D 1 [x„ - f(x;|c)]} 


AT A I 1 

— In |M| + - 5^[y„ -g(y„,x, l |b)] r M^ 1 [y„ -g(y„,x n |b)] (1.6) 

~ n— 1 

(L + M)IVln(27r/i). 


where x ra = x ”+l x " . Here Sd yn and S meas are the dynamical (first two terms) 
and measurement (next two terms) parts of the minus log- likelihood function. 
We note that Sdyn is the minus log-probability density in the space of dynam- 
ical paths and, in the limit of N — > oo, h — > 0, T = Nh = const, it coincides 
with the path-integral presentation obtained earlier in [21, 17]. 
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To find the general solution of the problem (1.1), (1.2) one can iterate 
optimization of S in the space of dynamical paths {x„} and in the space of 
parameters |c, b, D,m| (see [46]). From the view point of practical appli- 
cations for on-board IVHM system it is important to avoid heavy numerical 
computations and to restrict the solution to an analytical one whenever it 
is possible. It was shown in our earlier work that analytical solution can be 
obtained if prior PDFs for c, and b are chosen in the form of Gaussian distri- 
butions, and prior PDFs for D and M are uniform. 

Another key ingredient of analytical approach is successful factorization 
of the vector field. A general form of factorization that has proved to be 
very effective in many interdisciplinary applications can be written as follows 
[45, 46, 28, 30] , including specifically aerospace applications [37, 36] 


f (x|c) = F(x)c, g(y,x|b) = G(y,x)b, (1.7) 

where F(x) and G(y,x) are ((F x L) x L) and ((G x M) x M) matrices 
consisting of F (G) of diagonal blocks multiplied by 0*(x) (tpj(y, x)) base 
functions of factorization. The F- and G-dimensional sets of base functions 
{(/>,} and {tpi} are arbitrary and known. 

For the chosen above prior PDFs and factorized vector fields the analytical 
solution can be used to infer model parameters (cf with [46, 30]) 


<D> 


( c ) 

w*(D) 


3* (D) 


h N ~ l 

N X! [ x ™ ~ ^ n c 

n— 0 

H^(D)w*(D), 


N - 1 

h E 

n— 0 L 


F 1 D 1 x r 


N —1 

FnD-'F: 

n— 0 




(1.8) 

(1.9) 

( 1 . 10 ) 

( 1 . 11 ) 


where F n = F(x n ), and the components of the vector v(x) are 


L 

An(x) = ^ 
1=1 


( 9F; m (x) 
dxi 


m = 1, . . . , F. 


( 1 . 12 ) 


The parameters of the measurement model can be estimated using the con- 
ditions = 0 and d QM cas = 0, recovering the least square results in the 
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(M) 

1 N 

= N [ y " ~ 

77,= 1 

nr I T 

b y n - G n b , 

(1.13) 

(b) 

= ©^(M)z^(M), 

(1.14) 


N 



z x,y{ M) 

Yi — ]_ 


(1.15) 


N— 1 



e^y(M) 

= h ^ M " 1 G„, 

(1.16) 


n— 0 


where G n = G(y n ,x n ). 

Eqs. (1.8)-(1.16), coupled with the optimization procedure in the paths 
space, represent the general Bayesian framework for learning a nonlinear 
stochastic dynamical system from measurements that are corrupted by noise. 
Using this approach we can develop a method of fast on-line tracking of the 
time-varying parameters of non-stationary systems, as described below. 


1.3 The Lorenz system 

We start with the archetypical chaotic nonlinear system of Lorenz, 


ii = a (x 2 — Xi) + £i(t), | 

x 2 = r x\ - x 2 - xi x 3 + &(t) : > (1-17) 

X3 = X!X 2 - bx 3 + £ 3 (f), J 


augmented by zero-mean Gaussian noise processes £z(t) with covariance 
(£;(f) £i> (t')) = Dm 5(t — t'). Synthetic data (with no measurement noise) 
were generated by simulating (1.17) using the standard parameter set cr = 10, 
r = 28, b = |, and for various levels of dynamical noise intensities as ex- 
plained below. The phase portrait of the Lorenz system with dynamical noise 
is shown in Figure 1.3.1 along with the noiseless case to visually convey the 
difficulty of the inference problem. We now demonstrate that the formalism 
of dynamical inference outlined in the previous section allows one both to 
estimate parameters of the known dynamical equations and to to effectively 
perform model discovery when the vector-field of the Lorenz model is assumed 
to be unknown. We emphasize that both algorithms work well in the presence 
of strong dynamical noise. 



10 


Book title goes here 




FIGURE 1.1: The phase portrait of the chaotic nonlinear Lorenz system 
(1.17) with the standard parameters (see text): (a) deterministic system; 
(b) stochastic system with strong dynamical noise, simulated with a diagonal 
diffusion matrix having elements Du = 1500, D 2 2 = 1600, and D 33 = 1700. 
(All quantities in the equations and figures are dimensionless in this paper.) 


1.3.1 Parameter estimation with strong dynamical noise 

In parameter estimation, the functional form of the nonlinear force field, 
i.e. the right-hand side of (1.17), is assumed known, and the unknown values of 
the associated coefficients are then estimated from data. This is the approach 
reported in [19], where the diffusion matrix is taken in the form D = r 2 I, and 
the unknown parameters {a, r, 6, r 2 } are estimated via extensive numerical 
optimization of a cost function by simulated annealing and back-propagation 
techniques. We now demonstrate that our algorithm can estimate the param- 
eters of the system (1.17) analytically very efficiently and with high accuracy. 

First we notice that since the diffusion matrix is diagonal, our algorithm 
is reduces in this case to the trivial one-dimensional analytical solution of the 
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TABLE 1.1: Inference results for the 

parameters of the system (1.17) with strong 
dynamical noise. A synthetic data set of 4,000 
points was generated for each case by simulating 
the system with a diffusion matrix D = r 2 I, and 
subsequently sampling its trajectory with 
h = 0.002. 


Parameter 

Value 

Estimate 

a 

10.00 

9.9039 

r 

28.00 

28.3004 

b 

2.667 

2.8410 

T 

40.00 

39.9108 


problem for each equation in the form (cf with (1.8) - (1.11)), 

c, : = HA'w, : , i = 1 . . . 3 

where 


N-l , 
Wil — ^ ^ 

n — 0 ' 

T 2 d(j>u \ 

2 dxi ) 

N-l | 


■ 4 >il^iL 

11 

M 



n =0 ' 

\ ■ ■ 



Noise intensity is found according to (1.8). We note that in each equation 
we now have different basis functions (fu- For the first equation we have the 
following two basis functions: tpu = aq and <pi 2 = X 2 ■ For the second equation 
we have: foi = aq, f>22 = X2, and <f > 23 = X\X 3. And for the last equation we 
have: <j> 31 = aqaq, <p:n = x 3 . 

Thus there are a total of 8 unknown parameters to be estimated: a seven- 
dimensional coefficient vector c and the noise intensity r 2 . (Note that this is 
already more ambitious than what was done in [19], since we are attempting 
to estimate all model coefficients, including those that are equal to ±1.) 

The convergence of our scheme is so rapid that it is feasible to use the 
algorithm in real time on “streaming” data. To make a fair comparison we 
use the same number of data points as in [19]. As an indication of the inference 
accuracy, we quote in Table 1.1 results for data simulated with the standard 
Lorenz parameter set and two values of dynamical noise intensity for weak 
and strong cases. 

1.3.2 Model reconstruction with strong dynamical noise 

We now assume that the analytical form of the nonlinear force field of 
the system that generated trajectory shown in Figure 1.3.1(b) is not known 
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TABLE 1.2: Inference results for a representative set of 

parameters of the model (1.18), obtained using 200 blocks of 
600,000 data points each, sampled at ft. = 0.005. True and inferred 
parameter values are shown along with the corresponding error 
(relative and absolute errors for the nonzero and zero parameters, 
respectively). The inference error is below 1% for all parameters, 
and much less for most. 


Parameter 

Value 

Estimate 

Error 

an 

-10.0000 

—9.9984 

0.0161 

021 

28.0000 

28.0139 

0.0496 

«31 

0.0 

-0.0052 

-0.5180 

021 

10.0000 

9.9982 

0.0178 

O 22 

- 1.0000 

-1.0051 

0.5120 

023 

0.0 

0.0031 

0.3072 

033 

-2.6667 

-2.6661 

0.0196 

bill 

0.0 

0.0002 

0.0179 

&211 

0.0 

0.0002 

0.0238 

fr311 

0.0 

-0.0004 

-0.0401 

frll3 

0.0 

-0.0001 

-0.0111 

fr213 

- 1.0000 

-1.0004 

0.0446 

Dll 

0.2867 

0.2865 

0.0587 

Diz = D 3 1 

0.1069 

0.1061 

0.7657 


a priori. In this setting, it is more appropriate to refer to the inference prob- 
lem as model reconstruction. In practical terms, the main difference between 
parameter estimation and model reconstruction is in the number of unknown 
parameters involved, which is typically an order of magnitude larger in the 
latter case. To solve this problem within the formalism of dynamical inference 
one may adopt a parametric model of the form 

3 3 

xi = aw xy ( t ) + y: bwi"Xi*(t)xi»(t) + €i(t), (1.18) 

i'= i i',i"= i 

1,1' ,1" = 1,2,3. Including the elements of the (symmetric) diffusion ma- 
trix D, we now have a total of 33 unknown parameters comprising the set 
M. = {{a;/'}, {bwi"}, {Dui}}. Despite the restriction to linear, bilinear, and 
quadratic polynomial basis functions, (1.18) still represents an extremely 
broad class of dynamical models. Assuming no measurement noise for sim- 
plicity, the application of our algorithm entails the use of equations (1.9) - 

d-11)- 

The accuracy of the reconstruction depends on a number of factors. We 
have observed that it is generally possible to achieve arbitrarily accurate in- 
ference results with a (sufficiently small) fixed sampling interval by increasing 
the total duration of observation; this is true even in the case of a full (i.e., 
non-diagonal) diffusion matrix. Indeed, we were able to achieve highly accu- 
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rate parameter estimates for sampling intervals ranging from 1CP 6 to 0.01 and 
noise intensities ranging from 0 to 10 2 . As an example, we summarize in Ta- 
ble 1.2 our inference results for the model (1.18) with a full diffusion matrix, 
showing extremely high accuracy. 

Finally, we would like to demonstrate the importance of the Jacobian pref- 
actor included in our likelihood function by examining the inference results 
obtained with and without this term. Our analysis shows that the omission of 
the prefactor in the likelihood function results in a systematic underestima- 
tion of this parameter, whereas the inclusion of this term leads to an accurate 
inference as it optimally compensates for the effects of dynamical noise. 

The computational efficiency of our algorithm also allows us to lift the 
practical limitation on the total number of data points used for inference 
in previous work and to extend substantially the dimensionality of the model 
space. As a consequence it can be efficiently applied to deal with a more general 
problem of model reconstruction, when the functional form of a nonlinear 
vector field is unknown 

We now consider a selected set of examples of dynamical inference in 
aerospace applications. 


1.4 Three tank problem 

In section we describe an application of the dynamical system algorithm 
to a three tank [49, 27] which is formed by a sequence of three interconnected 
tanks and has been declared as a benchmark for fault detection and diagnosis 
in dynamical system [50]. The system schematic is shown in Figure 1.4. We 
consider this system in the context of ground support system for liquid fuel 
filling system. The key feature of the analysis presented below is that the 
underlying dynamics is hybrid system, i.e. the system exhibits both continuous 
flow and discrete jumps. Accordingly a of hybrid probabilistic approach is 
required for its inference. It is shown that the method of dynamical inference 
introduced above can accommodate abrupt changes in the model parameters 
and can be used to develop hybrid probabilistic algorithm. 

Accordingly the model equations are slightly modified to have the form 

hi = £isin(wt) - C 12 (£ [hi - H 12 \ - f[h 2 - # 12 ]) 

— A 12 \J fo [hi — Hi 2 ] — ci 3 (£ [hi — H 13 ] — £ [h 3 — H 13 ]) 

—Li3\/ fo [hi — H 13 } + f\ [^-l] + 0 ii£i (t), (1-19) 

hi = ~£ 2 sin (ut) + C 12 (£ [hi - i/ 12 ] -f[h 2 ~ H r2 }) 

+ C 23 (£ [h 3 — H 23 ] — i [h 2 — H 23 }) — L 23 \/ fo [h 2 — H 23 ] + 

h 3 = e 3 sin (ut) + c 23 (£ [h 3 - H 23 \ - £ [h 2 - H 23 ]) + cr 33 £i (t), 
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FIGURE 1.2: Schematic of Three-Tank bench Mark System. 


where 


£(f) = t ■ 0 (t); f 0 (t) = (1 + tanh(a • t))/ 2; 

h(t) = 0.2 • / 0 (0.2 — t) + (/ 0 (0.8 — x) — / o (0.2 — t)) + (1.20) 

0.2 • (/ 0 (1 - t) - /o(0.8 - t)) 


and 0(t) is the unit step function. Here H, t j represent liquid levels in each 
tank with respect to the connecting pipes, fi(t) is the known input flow, 
Cij represent parameters of coupling between the tanks, terms proportional 
to £i correspond to the vibrations of the tanks (in the present version only 
tank 1 is vibrating), terms proportional to L i3 correspond to the leaks, finally 
terms proportional to a 3 j correspond to the random vibrations of the tanks. 
Here £,(f) are white Gaussian noises and (J l3 are the amplitudes of random 
vibrations. 

The known input flow and the fi(t) and measured signals h\{t), /12(f), 
and h 3(f) are shown in Figure 1.4. The values of the known relative heights 
in these simulations were Hi 2 = 0.5, H 23 = 0.2, Hi 3 = 0.6. We emphasize 
that the time instants of the jumps in the flow signals h 2 (t) and /13(f) is 
not know and have to be inferred alone with other model parameters. The 
coefficients corresponding to the jumps are discrete parameters of our hybrid 
model. Note also that all three measured flow signals are perturbed by the 
periodic vibrations with known frequencies, but unknown amplitudes, and by 
random force. To apply algorithm introduced in ?? we have to specify the 
set of unknown variables 


c — {£1, £2, £3, C12, C13, C23, L \ 2 , L 13, L 23 , D 11, D 22 , D 33 } . 


(1.21) 
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FIGURE 1.3: Measured flow signals for the three tank benchmark system 
are shown by red solid line for hi(t), black dashed line for /i 2 (/), and by red 
dotted line for / 13 (f). The measured signals are compared with the input flow 
(black thick line) in the inset. The location of jumps for flows / 12 (f) and / 13 (f) 
are indicated by the arrows. 


and F-dimensional (F= 9) set of the base functions 

{(t>i{x, f)} = jsin(wt), (£ [hi - H 12 \ - £ [h 2 - H 12 ]) , \/£[hi - H 12 ], 

(£ [hi — Hif\ — f [h 3 — i/13]) , \/£ [hi — H13 ], (1.22) 

(£ [hi — Hi 2 \ — £ [h 2 — H 12 ]) , (£ [h 3 — H 23 ] — f [h 2 — H 23 \) , 

\f f [h 2 — H 23 ], (^ [h 3 - H 23 ] — f [h 2 — H 23 \ ) j" . 

It is assumed for simplicity that the measurement noise is small, which is 
justifiedin the case of application to the ground fuel filling system. Under these 
assumptions the inference of the unknown parameters of the model (1.19) from 
the measured time-series of the flow variables is a matter of straightforward 
extension of equations (1.8) - (1.11). In particular one has to check on each 
step of iterations if the base functions corresponding to the discrete parameters 
are identically zero and exclude related columns and rows from matrix H^(D) 
to avoid singularities in Eq. (1.9) 

An example of convergence of one of the model parameters is shown in 
Figure 1.4 (b). The circles show dynamics of the mean inferred value of the 
parameter C 12 . The bars show the dynamics of the corresponding standard 
deviation. The horizontal red line indicates the actual value of the model pa- 
rameter. It can be seen from the figure that the inferred value of C 12 approaches 
its true value as the total time of inference is increasing, simultaneously its 
standard deviation is decreasing indicating the convergence of the algorithm. 
Note also that initial value of coupling between tank 1 and 2 is inferred as 
zero. This value of the coupling parameter C 12 corresponds to the fact that 




FIGURE 1.4: (a) Number of operations as function of the total time interval 
of inference, (b) Convergence of the dynamical inference for the c \2 parameter. 
Time is shown in logarithmic scale. The horizontal red line shows the true 
value of the parameter ci 2 - 


at the initial time the second tank is closed and there is no flow between the 
two tanks. Next we can see that a jump to the correct value of C 12 = 2 occurs 
when the second tank begins to fill. Additional jump occurs when the third 
tank is also open (at around 2 sec). At this moment we can also observe a 
step-wise increase of the standard deviation indicating that the dimensionality 
of the system was increased in a step-wise manner. The number of operations 
required for convergence as a function of time is shown in Figure 1.4 (a) and 
is summarized in the Table 1.3. 

It can be seen from the table that discrete jumps in the coupling cij and 
leak Lij parameters are also detected and their continuous values are inferred. 
However, the convergence of the parameters is highly non-uniform. 

To clarify this issue let us make a few remarks. It was mentioned earlier 
that each block of data can be measured independently and used at the next 
step of inference provided that the results at previous steps are taken into 
account in the form of a prior distribution. In this case Eqs. (l.lO)-(l.ll) can 
be then written in the form (see [46]) 

Wfc = + hJ2 D -1 X„ - y- 

n&N k 

= 3fc - 1 + h F^D 1 F n - 

nGiVfc 

Applying this equations recursively one can show that the covariance matrix 
E*, of the posterior distribution is a sum over all the blocks and has the 
structure of a Kronecker product 

E fe = h<& ® D -1 , 


(1.23) 

(1.24) 


(1.25) 
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TABLE 1.3: Inference results for a set of parameters of the 

model (1.19), obtained using 300 blocks of 6000 data points each, 
sampled at h = 0.002. True and inferred parameter values are 
shown along with the corresponding standard deviation (STD). 


Parameter 

Value 

£l 

0.1 

Cl2 

2.0 

Cl3 

2.0 

L\2 

0.07 

L 13 

0.04 

L 23 

0.08 

C 23 

1.0 

D 11 

9.0 x 10" 

D 22 

4.0 x 10~' 


Estimate 

STD 

0.0993 

0.01 

1.988 

0.18 

1.999 

0.01 

0.047 

0.19 

0.068 

0.23 

0.08 

0.001 

1.004 

0.006 

8.9 x 10" 4 

1.4" 5 

4.0 x 10" 6 

1.3 x 10 


where 


( Vh,nt/h,n • • • 

. ■, ; 

■ ■ ■ VtB.nVtB.n 

Accordingly, all non-zero elements of this matrix grow in time as T = hN. 
On the other hand the second term in (1.24) remains finite for a fixed number 
of points in one block Nk- Therefore, a^-i approaches E/,. for large T and Ck 
becomes 



c fc « c k -i + D ® ^ fc 1 ^2 'xu-'y ■ 

n£Nk 


(1.26) 


To analyse the leading order contribution to the residuals of {ck} given 
by the second term in (1.26) we can assume that noise matrix is known and 
constant, the contribution from the sum can always be finite and small for 
small and finite Nk- The only factor that represent the accumulating effect 
of convergence is 4>). 1 . It is clear now that convergence of the residuals is 
inversely proportional to the sum of eigenvalues of <&k- Therefore the presence 
of large eigenvalues of 1 slows down the convergence of all coefficients {q,-}. 

It is also clear that the choice of the base functions has a strong effect on 
convergence the dynamics of the system. In particular, for polynomial base 
functions and small noise intensities the smallest elements of correspond to 
the highest powers of polynomials and usually correspond to large eigenvalues 

°l' */,’ 

Therefore, to achieve the best results in dynamical inference on-board it is 
important to identify the key non-stationary parameters for each fault mode 
and too learn most of the stationary parameters in a preliminary analysis 
of the system. Next, by incorporating real-time inference into the inferential 
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learning framework and excluding all but the most important non-stationary 
parameters from the tracking procedure one can improve the time resolution 
of the method by orders of magnitude. 

We next consider an application of the dynamical inference in the context 
of developing on-board FD&P system for solid rocket motors. 


1.5 In-flight decision support for SRMs 

Safe exploration of the space require development of the new algorithms 
with the overarching goal of extracting an information from available data 
with no “misses” and no “false alarms” . The existing algorithms rely heav- 
ily on the linearisation of the equations of motion and setting conservative 
margins for the nominal values of the control parameters. Yet a result the 
prognostic reliability of such algorithms is limited because of limited num- 
ber of sensors available on-board, limited thrust vector control authority, and 
short time window between the detectable onset of a catastrophic failure. It 
is very desirable to embed standard algorithms into a Bayesian framework to 
facilitate FD&P system. 

Here we demonstrate an application of the dynamical inference algorithm 
to a few fault scenarios including: (i) step-wise nozzle blocking with progres- 
sive steady burn-out, (ii) nozzle blocking with time varying fault parameter 
and neutral thrust curve modelling a possible “miss” situation, where sud- 
den increase of the pressure follows a prolonged period of small deviations of 
the pressure from the nominal value; (iii) and finally we will model a possi- 
ble “false alarm situation where pressure following a sudden increase in the 
chamber pressure, induced by a cloud of solid particles, returns to its nomi- 
nal value. To simplify the discussion we restrict our analysis to the case of a 
subscale motor. An extension of this approach to a large segmented motors 
can be found in [31]. 

1.5.1 Internal ballistics of SRMs 

The internal ballistics of the SRMs in the presence of the fault can be 
described by the following set of stochastic partial differential equations rep- 
resenting conservation laws for mass momentum, and energy of the gas (see 
e.g. [3, 39, 40, 31]) 

f dt ( Ap ) = -d x ( Apu ) + p p r h l + d^t), 

< d t {Apu) = ~d x ( Apu 2 ) - Ad x p + p p r b lu s + d 2 &{t), (1-27) 

[ d t ( Ape t ) = -d x {Aph t u) + Hp p r b l + d 3 ^ 2 {t). 

Here A is the port area, e t = cyT + u 2 / 2 and h t . = cpT + u 2 / 2 are the 
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total energy and total enthalpy of the gas flowing with velocity u, temper- 
ature T, pressure p, and density p. The propellant properties are described 
by the heat of combustion H , density p p and the burning rate given by the 
following expression r b = r c (p/p c ) n . To model performance of the functionally 
graded propellants we introduce the port perimeter l(t, x) as a given by design 
function of the burned distance R(x), e.g. for a progressive burning l = 2nR 
and for a neutral thrust curve l = const. The burning surface is given by the 
integral A b (t) = l(t,x)dx and, therefore, is determined by the propellant 
grading function. The increase of the volume of the combustion chamber dur- 
ing the time interval h is AV = A^^tfAR, where A R is given by the burning 
rate 

^ = r b (p 0 )p n + d 4 f 3 (t). (1.28) 

To model various uncontrollable sources of noise (such as cracks and case 
vibrations) that may become essential in off-nominal conditions a random 
component in the propellant density is introduced that results in the random 
forcing in Eqs. (1.27) and (1.28) with amplitudes di. Various fault modes in 
SRMs can now be modelled within the set of Eqs. (1.27) by choosing the time 
scale and direction of the geometrical alternations of the grain and case. For 
example, the bore choking fault can be modelled by introducing fault induced 
changes to the local port area A(x ); the crack dynamics can be modelled by 
introducing crack induced changes to an effective port perimeter; the nozzle 
blocking and the case breach faults can be introduced as changes to the nozzle 
throat area A t (see e.g. [32] and discussion below). 

1.5.2 Low-dimensional performance model 

To be able employ dynamical inference algorithm on-board one has to re- 
duce model (1.27) to a set of ordinary differential equations allowing for an 
analytical estimation of the model parameters and their fault induced devia- 
tions. To this end we introduce a low-dimensional performance model of SRM 
in the “filling volume” approximation [39, 40]. For a subscale motor the ratio 
of the gas velocity to the stagnation speed of sound Co is small {u 2 /cf'l) ev- 
erywhere along the propellant grain. Under these conditions one can integrate 
Eqs. (1.27) along the grain axis to obtain 

( d t {p 0 V) = - (puA) | L + p p r b (p 0 )A b + di£i (t),_ 

< d t (poeoU) = - (puAh t )\ L + Hp p r b (p 0 )A b + d 2 6(t), (1-29) 

{ d t R = n (p 0 ) , 


where subindex 0 corresponds to stagnation values of the flow parameters. 
Finally, taking into account the sonic condition at the nozzle throat in the 
form 


(Apu) \ L = PtUtAt , 


(1.30) 
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and introducing dimensionless variables (normalized by arbitrary reference 
values of pressure and density p s and p s ) 

P-^Po/Ps', P * Pol Ps (1-31) 


we arrive at the following from of the low-dimensional performance model 
(LDPM) of the SRM operation in nominal and off-nominal regimes 


dp 

dt 

dp 

dt 

dR 

dt 


c 0 "fTA t (t) 

V(t) 

c 0 TA t (t) 


pVpTp 


A b (t)r b 


Vpp- 


V(t) 
A b (t)r b 


(7 P P ~p)p n + 


d-2 (7 ~ 1) 

V(t) 


V{t) V(t) 

r b p n + d 3 (,3 (t) , V = A b (t) r b p n . 


( Pp - P) P n + 


V(t) 


Ut), 

(1.32) 


Here the burning area A b {t) = f(R(t)) is given by the design of the propellant 
grading and is a function of R. It is now possible to characterize the dynam- 
ics of the SMB operation by the parameters averaged over the length of the 
combustion camera. For example, the cracking and bore chocking faults will 
affect A b (t) and V(t), while the nozzle failure or the case breach faults may 
affect A t (t), the growing intensity of the noise can also be an important indi- 
cation of the fault. The dynamics of the fault in these settings is characterized 
by an unknown time-dependent fault function ff au lt(t) that in many cases of 
interest can be assumed to be a polynomial function of time with unknown 
coefficients. We now demonstrate that abrupt changes of the model parame- 
ters can be accommodated within introduced Bayesian inferential framework 
of on-board FD&P system for SRMs. 


1.5.3 Estimation of the parameters of nozzle blocking 

In our first example we consider inference of the nozzle blocking fault 
modelled as a step-wise change in the nozzle throat area A t . An additional 
rationale behind this example is that the significant step-wise change in the 
parameters of nonlinear models is one of the most plausible reasons for the 
failure of alternative Kalman-fflter based algorithms. The nozzle blocking is 
modelled as a sudden change of the nozzle area to the effective value A et = 
A t + A At (cf. [32]). The corresponding geometry change of the nozzle is shown 
in Figure 1.5.3(a). To estimate the change of the nozzle area we notice that 
it is the coefficient a = c ° 7 yp) ? which is proportional to A t . We also notice 
that time-series R(t), V(t), and A b (t) are completely determined by the time- 
trace of pressure p{t) and can be excluded from the inference algorithm. The 
results of dynamical inference obtained using this algorithm are shown in 
Figure 1.5.3(b) and summarized in the Table 1.4 

It can be seen that the dynamical inference algorithm provides an accurate 
estimation of abruptly changing SRM parameters and that parameter’ PDFs 
become sharply peaked about the true value. 
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TABLE 1.4: The results of the parameter estimation of the model (1-27). 

The total time of the measurements in this test was f=0.5 sec, the sampling 
time was /i=0.001 sec, and the number of measured points was A=500. The 
propellant grain has cylindrical geometry. 


Parameter 

Value 

Estimate 

Relative Error 

27 P P 

2575.34 

2572 

1.1% 

-{coiTA t )/V 

-1101.75 

-1103 

1.1% 

d\ 

0.0902 

0.0906 

0.4% 

4 

0.0902 

0.0906 

0.4% 




FIGURE 1.5: (a) Geometry of the nozzle blocking model. The contour of 
the combustion chamber before the fault is shown by the black solid line and 
after the fault by the red line. The grain location is shown by the dashed 
line, (b) Estimation of the value of the parameter a before (left curve) and 
after (right curve) the fault. The dashed line shows the actual value of the 
parameter. The solid lines show the PDF of the parameter estimation with 
t= 0.5 sec, /i=0.001 sec, iV=500. 


1.5.4 Predicting “misses” in the fault detection 

To illustrate prediction of “misses” in on-board FD&P system let us con- 
sider a situation when small pressure induced deviation from the nominal 
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value persists for a few second prior to the rapid crossing the “alarm” level. 
Such fault dynamics is typical for combustion instabilities in solid propellant 
motors [12] and represents an archetypal example of “misses” . 

This situation is illustrated in Figure 1.5.4, where measured pressure signal 
(black solid line) crosses the alarm level (dashed line) and this initiating the 
alarm at approximately tA = 15 sec. The overpressure fault occurs at ip = 
17 sec and the time window between the alarm and a “catastrophic” event 
becomes too short, which can be considered as a model of “miss” situation. 
We model the “misses” by introducing fault-induced changes to the nozzle 
throat area and assume that the time evolution of the nozzle fault is highly 
nonlinear and can be described by a polynomial function 

A e t = A t0 — (ar + /3t 2 + 5r 3 ) . (1.33) 

This polynomial approximation of the fault dynamics is substituted into Eqs. 
(1.32). It is the presence of high-powers of time the cause rapid deviation 
of pressure from the nominal value at tp = 17 sec. The earlier detection of 
“misses” in this context is reduced to a two-step procedure: (i) detection of the 
fault initiation and (ii) earlier detection of high-power terms in the polynomial 
fit to the time evolution of the fault. Accordingly, the pressure time-traces have 
to be followed continuously on-line and fast algorithm of dynamical inference 
has to be employed for on-board estimation of the fault parameters. 

To simplify further discussion we neglect random term in the equation for 
the burn distance and assume that the thrust curve is neutral. Then burn 
distance, burning area, and the combustion chamber volume are a known 



FIGURE 1.6: Example of possible time variation of the pressure fault (black 
line) representing a possible “miss” situation. The blue dashed and red solid 
lines indicate the “alarm” and the “catastrophe” levels respectively. Note that 
the time window between the “alarm” and the “catastrophe” is too short. 
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p, atm 


FIGURE 1.7: (top) Time evolution of the pressure build up after the nozzle 
blocking fault is shown by the back solid line. The results of the predictions 
build 1 sec and 2.1 sec after the fault are shown by green and blue lines 
correspondingly. Examples of the pressure time-traces predicted ahead of time 
using inferred values of the fault parameters are shown by the jiggling lines. 
The green lines bounding yellow shaded area indicate standard deviations of 
the predictions after 1 sec of inference. The blue lines bounding blue shaded 
area indicate standard deviations of the predictions after 2.1 sec of inference. 
The time moments of the predicted overpressure faults used to build the PDF 
of the case burst times as shown by the black circles on the red margin line. 
Fault occurs at t = 9 sec. (bottom) The PDF of the predicted values of pressure 
at t=14 sec build 1 sec (green line) and 2.1 sec (blue lines) after the fault. 
The dashed vertical line shows the dangerous level of the pressure. 


functions of the pressure time-trace 


R(t) = rb f p n {t')dt'\ Ab(t) = Abo = const; 
J o 


V{t) = V 0 + A b0 R{t) (1.34) 
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and the equations of the SRM internal ballistics take the form 


dp 

dt 

dp 

dt 


~ + ~ ^ + ( L35 ) 


where we have introduced dimensionless time t — > Tbt, appropriately scaled 
constant unknown noise intensities d\ and d' 2 , effective throat area A et (t) given 
by Eq. (1.33), and R'(t) = Vo/Abo + R(t). 

To apply an algorithm of dynamical inference we rewrite the factorized 
vector field in the form f(x,t) = C<f>. The set of the base function can be 
chosen in the form 


im} = % 



P P R 1 > VpPr' i y/PpR. > VpPr' } • 


And the matrix C can be written as follows 
C = 

"fp p —1 0 —ay — aya' —ay/3' — ajS' 0 0 0 0 

p p 0 — 1 0 0 0 0 —a — aa! —a/3' —a5' 

Here a = (c$T A t $) / (rbVo) and coefficients a', [3\ and S' are equal to the 
coefficients a, f3, and S in Eq. (1.33) scaled by Abo/Vo. 

The results of the detecting the “misses” are shown in Figure 1.5.4. To de- 
tect fault initiation parameters of the system are monitored in real time. Once 
fault is detected at time t,i the algorithm is continuously updating the inferred 
values of parameters estimated on increasing intervals of time elapsed from 
td . These values are used to generate a set of trajectories predicting pressure 
dynamics and to calculate PDFs of pressure ahead of time and of the instants 
of time corresponding to the case burst. Example of trajectories generated 
for two different time intervals 1 sec and 2.1 sec are shown in Figure 1.5.4 
(a) by green and blue lines respectively. Examples of the PDFs of the time 
moment of the overpressure fault are shown in Figure 1.5.4 (b). It can be 
seen from the figure that predicted PDFs converges to the correct value 2.1 
sec after the fault extending the time window between the “alarm” and the 
fault almost three folds (cf. with Figure 1.5.4) and reducing the probability 
of the “misses”. Further detailed discussion and examples of application to 
an analysis of the results of the ground bring test and to a case breach fault 
diagnostic and prognostic in a large segmented SRM are give [29] . 
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1.6 Diagnostics of space vehicle stage separation failure 

The stage separation failure may have various origins including e.g. failure 
of the accelerating/decelerating motors during the separation. The combined 
expected failure rate is relatively high and is the second most common cause 
of launch failures [6] . Despite this fact a reliable in-flight diagnostic and prog- 
nostic system for stage separation failure is currently unavailable. Below we 
report a progress in development of such system. Development of an in-flight 
FD&P system for the stage separation failure is a challenging engineering 
problem [41]. The main difficulties stem from the fact that the phenomena 
underlying fault dynamics are highly non-linear and transient in structural, 
thermal, and fluid dynamical domains. Accordingly, the solution requires 3D 
analysis of the thermal/fluid/structure interaction in the supersonic flow. In 
addition the number of sensors available on-board is severely limited and the 
safe time window between the detectable onset of the fault and possible catas- 
trophic failure is typically a few seconds. To overcome these difficulties we use 
a combination of the physics based analytical model and high-fidelity simula- 
tions using ABAQUS [1] and FLUENT [2]. 

We investigate a problem of the second stage engine nozzle damage due to 
the impact with the first stage of the vehicle in the process of stage separa- 
tion. We describe high-fidelity method for analytical estimations of structural 
dynamics of the second stage engine nozzle, analyze nozzle damage due to the 
re-contact of the stages and discuss the results of the computer simulations of 
re-contact process. We describe method of nozzle damage diagnostics in real 
time using vehicle sensor data. 


1.6.1 Nozzle extension impact dynamics 

We consider two-stage, vertically stacked space vehicle. The vehicle is pro- 
pelled by the first stage booster at the first phase of the flight. After the rocket 
booster completes its mission, the stages separate and the second stage engine 
provides the propulsion during the second phase of the flight. Due to separa- 
tion faults, the first stage can collide with the second stage engine nozzle and 
damage it. As the first task we perform the analytical description of nozzle 
damage after the impact using finite element model and the theory of shells. 
Dynamics of the second stage engine nozzle extension under the impact is well 
approximated by Donnell’s shallow shell theory. Dynamics of the second stage 
engine nozzle extension under the impact is well approximated by Donnell’s 
equation of shallow shell theory [26] 


DA 2 w-V 2 r F 


, dw 

ch ln 


, d 2 w 
ph -» = 


f(t,s,e ) 


(1.36) 
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^A 2 f + V> = 0, (1.37) 

Gr 

where w transverse displacement, F - an Airy type of stress function 

_ ld_ 1 d 2 2 _ 1 d 2 

ds 2 s ds s 2 sin 2 a dO 2 ’ R s tan a ds 2 ’ 

D and G are bending and shear stiffness c - damping coefficient, h - thickness 
of the plate. General view of the cone geometry is presented on Fig. 

Consider the small end is the clamped one and the large end is free. As- 
suming this case, the boundary conditions are 


diu 

u = v = w = — — = 0, at s = Si 
os 

N s = s s g = V s = M s =0, at s = s 2 - 


(1.38) 

(1.39) 


u, v, and w are the orthogonal components of displacement in the s, 9 , and 
normal directions, respectively, N s ,S s g are the membrane forces related to 
the Airy stress function, V s is the Kelvin-Kirchhoff shear, M s is the meridional 
moment resultant [26]. 

Let us show how eigenfunctions and eigenfrequencies can be found. Con- 
sider separation of variables 


w(t,s,6 ) = e lbJt w(s) sin(n0), 

F(t,s,6) = e lult F(s) sm(n9). 

We obtain for the system (1.43) (1.43) 

/ Da(n) —b \ / w(s)\ _ ( hpu; 2 \ / w(s)\ 

^ b G~ l a(n) ) \f{s)) ~ V 0 ) \ F { S )) 


(1.40) 

(1.41) 


(1.42) 


where a = 


d t 

ds 4 


2 d 3 
s ds 3 


l+2n 2 d 2 , l+2n 2 d 

s ds 2 ' s 2 ds 


n 2 (4— n 2 ) ] 


,6 = 


1 d z 
s tan a. ds 2 ' 


The eigenfunctions w mn (s) , F mn (s) satisfying (1.42), (1.38), (1.39) and nat- 
ural frequencies ui mn of the cone can be find by Galerkin’s or Rayleigh-Ritz 
methods. Expanding the radial displacement w, and Airy function F(t,s,0 ) 
in the series of the shallow shell eignemodes 


w(t, s, 9) 
F(t, s, 9) 


y ^Cmn(t)w mri (s)sin(nQ), 

mn 

Y. dmn ( t)F mn (s)sin(n9 ) 

mn 


(1.43) 

(1.44) 


Using that normal modes with different indexes are orthogonal we obtain 
equations governing the amplitudes c mn as 
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d 2 

Cmn 0 r dc 

mr 

where 5 = c/2p and o; nm are frequencies of the eigenmodes. 


“1“ mn)C-mn fmn (^) 


(1.45) 


fmn{t) = j^ j j f{t,s,6)w mn {s)sin{n6)dsde 


To analyze response of the nozzle to the impact, we build finite element 
model of the nozzle extension using ABAQUS software package. Because the 
nozzle shell is much thinner than the impacting part, the later is modeled as a 
rigid body. We fix nozzle extension at the base and use dynamic explicit mode 
with general contact properties to simulate the impact. Typical results of the 
simulations are shown in Fig. 1.6.1(a). It can be seen from the figure that the 
nozzle is damaged during the impact. The damage can be characterized by 
two key factors: the area S of the bended part of the nozzle, and the angle of 
the bending. Next the geometry of the nozzle damage obtained in ABAQUS 
simulations is translated into the FLUENT model of the flow through the 
damaged nozzle and the torque amplitude induced by the flow on the damaged 
nozzle is estimated. As a result of this analysis we can use the inference of the 
impact strength to predict damage size and the torque and thrust imposed by 
the flow on the damaged nozzle. 

To analyze in more details the nozzle damage dynamics we use the solution 
of Donnells equations Eq. 1.36 and Eq. 1.37 to simulate the impact-induced 
torque applied to the nozzle and perform computer simulation of the upper 
stage thrust vector control response to the impact. 



FIGURE 1.8: (a) Nozzle extension damage during stage separation failure 
obtained by simulations of the high-fidelity model of the nozzle extension 
impact in ABAQUS. (b) Predictions of the flow distribution in the damaged 
nozzle using FLUENT simulations. 
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FIGURE 1.9: (a) Nozzle angular positions tilt (red) and rock (blue) during 
the impact, (b) Induced by the impact actuators forces components tilt (red) 
and rock (blue) as functions of time, (b) Induced by the impact actuators’ 
forces components tilt (red) and rock (blue) as functions of time. 

1.6.2 Second stage thrust vector control system operation 
simulation 

Second stage thrust vector control system consists of three integrated com- 
ponents: turbine pump assembly, hydraulic distribution system and actuator 
assembly. The gas from main propulsion system starts the turbine pump as- 
sembly. Turbine pump assembly produces power for hydraulic pump. The 
pump circulates the fluid through the hydraulic distribution system. Two hy- 
draulic supply strings support two actuators, which change the nozzle attitude. 
At nominal regime of stage separation the actuators do not operate and nozzle 
attitude does not change. When the first stage collides with the nozzle, the 
nozzle changes its angular orientation and the control system respond to the 
impact. 

For the analysis of nozzle damage due to stage separation we use dynamic 
model of the vehicle second stage thrust vector control. We simulate the impact 
torque for a given stage separation fault and implement the impact torque into 
the second stage thrust vector control system simulation software. In the result 
of the simulation we obtain nozzle attitude changes and the actuators’ forces 
as functions of time. The nozzle attitude changes during the impact are shown 
in Fig. 1.6.2(a) and the actuators’ forces are shown in Fig. 1.6.2 (b). 

We simulate impact torques for several separation faults, obtain impact 
torques for each fault and simulate thrust vector control response to each 
case. Different separation fault cases are described by different amplitudes of 
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FIGURE 1.10: Actuators tilt (red) and rock (blue) torques components 
obtained with the MATL AB/SIMULINK simulations of the actuator forces 
as functions of time. These components coincide with torques components 
obtained in the ABAQUS simulations. 


torque oscillations and, as the result, different nozzle angles and actuators 
forces response. 


1.6.3 Diagnostics of stage separation faults 

The objective of this section is to describe how we can detect separation 
faults and estimate nozzle damage using the second stage engine sensor data. 

1- st Method: Second stage nozzle attitude analysis 

We assume that we simulated several possible separation faults and obtained 
nozzle attitude changes for each case. If we have real time telemetry data of 
nozzle angular orientation in real time we can compare the telemetry data 
to each case from the previously simulated cases and to pick up the case of 
maximal correlation with the telemetry. Since we already studied the cases 
and estimated consequences of the impact, we can provide the prognostic of 
the investigated case. 

2- st Method: Impact torque analysis 

The second stage engine nozzle attitude changes are given by the equation 

J(p + Bip + Kip = Q (1-46) 

where J, B , I\ are the coefficients, <p is (rock/tilt) nozzle angle. The total 
torque applied to the nozzle Q = Q a + Q c is the sum of the actuator torque 
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Q a and the impact torque Q c . If we know the actuator force and the nozzle 
angle, we can compute the impact torque with the equation 

Q c = Jip + Bip + Kip — F a ■ l (1-47) 

where F a is the actuator force and l is the actuator lever arm. Measuring the 
nozzle attitude angles ip and actuator forces F a by the second stage engine 
sensors we can use the equation 1.47 and compute the impact torque applied 
to the nozzle extension in real time. The result of this simulation is shown in 
Fig. 1.6.3. 

The information about impact torque behavior will allow us to estimate 
the nozzle extension damage and to predict the consequences of the impact. 
We can compare the obtained impact torque to the result of simulation for 
possible separation fault scenarios and choose the most probable one based 
on the maximal correlation of the torque Fig. 1.6.3 with the simulated cases. 

It can be seen from the discussion above that the once the impact induced 
torque is detected the prediction of the damage is reduced to an analysis of 
the buckling eigenmodes of the thing conical shell that can be reduced to an 
analysis of a set of forced oscillators. 


1.7 Normal mode expansion method dynamics for local 
damage detection 

Metallic and composite plates are common members of airspace structures. 
As the size and complexity of space hardware grows, structural weight con- 
trol becomes crucial. Weight requirement usually involves the use of lighter 
weight structures such as composites. All composite structures have one basic 
handicap in common: they all share a failure mechanism that is so insidious 
as to make it difficult to impossible to discern when and where a serious flaw 
may occur. The reason is that structural failures will occur internally, out of 
view of normal visual means of inspection. Internal cle-bonds, de-laminations, 
cracks, and/or buckles are typical. For heavily loaded composites, such as the 
heavy lift launcher will undoubtedly employ, this type damage represents a 
major safety concern especially on man rated vehicles. The large area compos- 
ite structures like Payload Shroud, the Interstage, the Core Intertank, Storage 
Fuel Tank, Crew Composite Modulus etc. are under broad investigation to use 
them in Space Vehicles. Some of the structures peculiarities of these structures 
you can find in reports [9, 10, 44]. 

In this section we report some properties of the wave propagations in 
plates with application to development structural health monitoring (SHM) 
methods. Knowing the dynamics of the structures is important for maintain- 
ing structural integrity, safety, health and IHM as whole. In this section we 
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by using classical Kirchhoff thin plate vibration theory will study structural 
properties and show how these results can be used for structural health moni- 
toring. We will base our consideration on seminal work of Laissa (1969) where 
presented also references of the original papers. 


1.7.1 Governing equations 

Assuming the mid-surface of a plate as a reference, we shall use in-plane 
Cartesian coordinates. Let w be the deflection of the middle surface of the 
sandwich plate and than the governing differential equation of the isotropic 
thin plate can be presented from the classical thin plate theory. The trans- 
versely vibrating plates usually consider following equation [26] 

d 2 w 

DA 2 w + ph~g ^2 = f(t, x, y), (1.48) 

where D is the flexural rigidity (Bending stiffness), and p, h are density and 
thickness of a plate. Operator A 2 = V 2 V 2 - is biharmonic operator, f(t,x,y) 
term is lateral force per unit area responsible for external load or PZT excita- 
tion. In the case of thin homogeneous plate D = Eh 3 / 12(1 — v 2 ) is the flexural 
rigidity of the plate; E, v Young’s modulus, Poisson’s ratio, respectively. In 
the case of sandwich panel stiffens is [4] 


D = D S 


M , E c tj 

6 12 


Eftf(tf + t c ) 2 
2 


2 Df + D c + Dq. 


(1.49) 


Here Young modulus for facesheets is Ef and modulus of the core E c ,tf and 
t c are the thicknesses of the core and sheets, correspondingly. The sandwich 
plate has a thickness h = 2 tf + t c . Each term in the right hand side is denoted 
by its own capital letter D, p = (2pftf + p c t c )/h, pf,p c are the densities of the 
facesheets and the core. In expression (1.49) we have Df <C D 0 , D c <C D 0 . 

The displacement field w(x,y,t ) = Woexp(i(u>t— kx)) describes transverse 
plane wave traveling in the x direction. As a result, solving (1.48) for the 
wavenumber we get two real roots which describe two wave propagating in 
opposite directions: 

k = ± \/L0 2 ph/D, (1.50) 

while the two imaginary roots correspond to evanescent fields. The phase 
speed v p of the flexural wave with frequency w in this case is function of the 
plate stiffness and mass per unit area 


v p = aj/k = ± \J \o 2 D / ph. 


Let us consider generation and propagation of the flexural waves based 
on results of the eigenfunctions and eigenfrequencies obtained for rectangular 
plate. The collection of solutions for rectangular plates can be represented by 
deflection functions w(x,y) as the product of beam functions 
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w(x,y) = EE c mn X m (x)Y n (y), (1.51) 

m n 

where X m (a:)and Y n (y) are chosen as the fundamental mode shapes of beams 
having the boundary conditions of the plate. This choice of functions then 
exactly satisfies all boundary conditions for the plate, except in the case of the 
free edge, where the shear condition is approximately satisfied. The possible 
boundary conditions along the edges x,y = 0 and x = L x ,y = L y determine 
mode shapes X m (x),Y n (y). In this case the eigenfrequencies u> are given by 
the formula [26] 


W 2 = 


n 4 D 

Wp 


Gt + Gy-J + 2— | \vH X Hy + (1 - V)JyJy \ 


V T 4 

L v 


L l 


(1.52) 


where G, H , and J, are given functions of natural numbers determined by 
specific boundary conditions (Table 4.1 [26]). As it was noted by Leissa there 
are 21 combinations of simple boundary conditions (i.e. , either clamped (C), 
simply supported (SS), or free (F)) for rectangular plates. Frequency param- 
eters are expressed in terms of coL x \/phD, where L x is a length dimension. 
For example for SS boundary conditions for each edge we have 


w(x, y) = 0, at x = 0, L x , y = 0, L y . 


The problem with Simply Supported edges is one of the most simplest for 
analytical calculation and for the rectangular plate 


w{x,y) = EE' 


, sin 


m 'KX ^ ^ mry ^ 


V L x 


V Ly J 


(1.53) 


where L x ,L y are the rectangular plate dimensions 

G x = to, H x = J x = to 2 , G y = n, H y = J y = n 2 (1-54) 

and as a result, from (1.52) we have next eigenfrequencies 


^ mn 


tt 4 D 

ph 



where m, n = 0, 1, 2..., to = n yf 0. 


(1.55) 


1.7.2 Evolution under PZT actuator 

The evolution of the plate under external load can be find by considering 
amplitudes in equation (1.51) as a functions of t. Substituting (1.51) into 
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(1.48), multiplying both sides by X m i(x)Y n i(y) and integrating over the area 
of the plate we obtain 


^'^[{ C rnn( t ) 


L X Ly 


I ^ 

I Cmn^n 


X rn (x)Y n {y)X rn '(x)Y n >{y)dxdy 


o o 


L x Ly 




j J f{t,x,y)X m '(x)Y n '(y)dxdy. 


(1.56) 


o o 


Because the modes are orthogonal left hand side is not equal to zero only for 
n = n',m = m! . As a result, system (1.56) can be presented as system of ODE 
for amplitudes c mn 


Cmnif) 4” ^mn^mnit) — fmnif ) 


(1.57) 


where 


Ly 


fmnif) = J J f(t,x,y)X m {x)Y n (y)dxdy. 


0 0 


(1.58) 


We can easily solve the system (1.57) analytically when surface forcing is 
described by Dirac delta function <5(.) in space and harmonic in time. In this 
case 

fit , x, y) = cosiutfSix - x s )Siy - y s ), 


where point ix s ,y s ) is the position of source term on the plate. Introducing 
notations 

fmn(t) = foit)fmn, (1-59) 

where term f mn = L 4 L X m ix s )Y n iy s ) determines localization of the source 
on the plate, and foft) = cos(uA) we can write down the solution 


Cmn 


it) 


Co COS iajjnnt + of) + fmn 


cos iait + (3) 
- w 2 ) ' 


Amplitude Co and phase a are determined form initial conditions. Let us con- 
sider c mn = Cmn = 0, than Cq = — f m n c ° s ^ and solution for amplitudes 
and phase are 


Cmn{f) — 


frr 


i^mn - <X 2 ) 


, .. cos /3 , . 

COS (cot + 13) COS (uJmnt + «) 

COS (X 


(1.60) 


a = tan 


-l 


LO 


LUrt 


tan(/3) 


(1.61) 


Solution evidently explains one of the possible ways of SHM. The denominator 
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is equal to zero at w = oj mn which corresponds to resonances of the plate. 
Any structural defects will change u) mn and by sweeping frequency u> we can 
determine shift of u) mn . 

Let us consider next example when pulse is presented by the form [25] 

fo = e~ t l T cos cot. (1.62) 


The excitation signal (1-62) makes it possible to present solution of the system 
as [25] 


Cmn 


frr 


m 2 mn cos u mn t + (TLUmn) + 2 u 2 ) sin gj mn t 

ro mn + 4w 7' r2 


f tlr w hn cos - 2 (w/r) sin ojt 

+Jmn '-4 w 2 /r 2 


VDZ, 


(1.63) 


where w 2 mn = u 2 mn + 1/r 2 - w 2 . 

Let us consider how wave propagating from point source can interact with 
local damage. 


1.7.3 Model of the damage 

Local stiffness reduction can be caused by different faults. The loss of 
bending stiffness can be caused by debond, delamination, honeycomb crash or 
reduction of the modulus of facesheet due to impact. Any of these structural 
failures will lead to change of the stiffness and mass parameters. As a result in 
framework of plate theory we can generalize damage as domain with different 
structural properties than rest of the plate. 

Let us consider damage model as a decrease in structural local stiffness by 
function g(x — Xd, y — yd ) which reflects changes in the local structural material 
elastic coefficients. Term in the left-hand side of equation (1.48) is modified 

d 2 w 

Dg(x -x dl y- yd)A 2 w + P h ^^ = f(t, x, y), (1.64) 

where function describing damage is centered g{x — Xd, y — yd) at the certain 
point (xd, yd) characterizing local change in stiffness. We transform equation 
(1.64) to the form 

d 2 w 

DA 2 w + ph- 7 ^ = sf(t, x, y) + De(x, y)A 2 w, (1.65) 

where right hand side value e(x, y) = 1 — g(x — Xd, y — yd) is not equal to zero 
in a small domain of the damage. 

Considering that the value e(x,y) in (1.65) is sufficiently small we can 
apply perturbation technique: we find the solution for e{x,y) = 0 and then 
substitute displacement in right hand side in the expansion through normal 
modes form (1.51). In this case damage term works like a source of secondary 
waves generated in the structure. 
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Applying the same procedure of the expansion solutions system of equa- 
tions for mode amplitudes is 


Cmn{t) Cmn(t')^mn fmnif} A ^A rnn{x d ^yd) • 


where for calculation of the term eA ?I 


p A — 


LxLy ^ 


J-'x 

Y. J J e(x,y)ul l c k iX rn (x)Y n (y)X m >(x)Y n >(y)dxdy (1.67) 


we have to know exact distribution of the value e(x,y) 

As a result, we obtain solution for transverse displacement w(x, y, t ) which 
can be presented as 


N x ,N y 

w{x,y,t) = E ( t ) + fA mn (x d ,y d ) [1 - cos(2wf)]} x 

m,n = 1 



( 1 . 68 ) 


where c mn (t) is determined by expression(1.63) and eA (x d ,y d ) mn . 

In order to calculate (1.67) we consider c k i(t) at the damaged region as a 
given from solution for e = 0. Let us consider simplest model of damage when 
in local square we have a reduced stiffness by value e = const 


f £ for x d - D < x < x d + D, x d - D < y < x d + D, . , 

\ 0 for the rest of the plate. 


In the case of simply supported boundary conditions we have sine eigen- 
functions (1.53) and damaged term in (1.66) can be expressed as 

E w fcZ c M 

k l 

£A mn(x d ,y d ) = 2 — — [a km (x d ) - a krn (y d )][ai n (x d ) - ai n (y d )\, (1.70) 

7 T 

where a km {x d ) = ^ sin(fc - m) a km (y d ) = ^ sin(fc + ?n)^|^+^, 

ain{x d ) = (^sin(Z-n)gl^, a ln (y d ) = ^ sin(Z + n)^|^+j. 

Solution is given by (1.68) and these results are plotted on figure 1.11 
a) and b) where first graph represent total solution and plot (b) reflects a 
difference in displacements for healthy structure and damaged one. 

For computer simulation the force excitation point S is kept fixed at the 
center of the plate {L x / 2, L y / 2) and the damage is located at certain distance 
(point x d ,y d ) from the corner (0,0). Some of the results for transverse dis- 
placement are plotted in figure 1.12 as a function of ( x,y ) for different value 
of t and the simply supported boundary conditions along the edges. In the 
simulation we used Hanning windowed signal with number of pulses N 




36 


Book title goes here 


t= 6.31 58e-005 no J ANALYTICAL 


t= 0.00014211 differnce between analitical 



0 0 


0 0 


(a) 


(b) 


FIGURE 1.11: Analytical solution for w out of plane displacement for: t = 
0.0002sec - (a), and difference between the damaged and pristine signal - 
(b)(xd = yd = 0.5 to, d = 0.025m, x s = y s = 1 to, L x = L y = 2m, N x = N y = 
20, Input signal is determined by formula (1.63), r = O.OOOOlsec, / = 10 kHz, 
e = 0.1, uj = 2nf. Structural parameters are Ef = 10 11 Pa, pj = 1580 kg/m 3 , 
p c = 110 kg/m 3 , E c = 2 • 10 8 Pa, t c = 0.03m, t = 0.003m. 


where cj is the central frequency and signal is comprised of component waves 
over a close range of this frequency. Changing of the base frequency u> changes 
the wavenumber (1.50) and wavelength of the excitation signal. The influence 
of a local change in stiffness (different stiffness for the size 5 cm) on the 
propagation of flexural waves is considered. The simulation was made for a 
different number of the modes N x , N y . The mode number affects the accuracy 
of calculation. We simulated the structure till reflection from the boundaries 
is negligible. By increasing time of simulation incident pulse interfere with 
reflected one and patterns of wave distribution are much more complicated. 
The local stiffness lost changes the amplitude of the scattered signal but does 
not change general picture. We can observe that damaged region plays the 
role of the scatterer generating scattering waves. 

For the example expressed by (1.68) two distributions are plotted just to 
show evolution of field w(x, y). Clear indication of damage was shown for the 
time when wave reaches the damaged region (Figure 1.12c). The difference 
in the w(x,y) between the damaged and pristine signal is depicted in Fig- 
ure 1.12d). (distribution w(x,y) without first term in (1.68)). This makes it 
possible to characterize the damage. 

This study clearly shows that normal mode expansion method is capable of 


fo(t ) 


f sin (jjt) sin (cot) , for t < to = N-jt/cj, 
\ 0, for t > t 0 - 
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FIGURE 1.12: Out of plane displacement w(x,y) at t = 2 • 10 _5 sec- (a), 
t = 3-10 ~ 4 sec- (b) for nondamaged plate, for damaged one at t = 2-10 ~ 4 sec- 
(c) and difference between the damaged and pristine value of w(xpy) at t = 
2 • 10 _4 sec- (d). Input signal is Hanning windowed signal with N = 3.5, 
/ = 10 kHz, Structural parameters are Ef = 10 n Pa, pf = 1580 kg/m 3 , p c = 
110 kg/m 3 ,E C = 2 • 10 8 Pa, t c = 0.03m, tf = 0.003m. 
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capturing the physics of wave propagation in plate structures. We see a scat- 
tering pattern which in general closely resembles that interaction of wave with 
damaged region. Using the results obtained earlier, it is analytically shown and 
verified by computer simulation that the normal function expansion method 
describes damage state of the plate by changing the evolutionary dynamics of 
several modes. As a result, dynamic inference method is very promising for 
SHM of the plates. Taking into account large number of eigenmode expan- 
sion increases the accuracy but at least ten terms is sufficient to grasp main 
characteristic features of the damage presented in the plate. 


1.8 Dynamical inference of a set of coupled oscillators 

It was shown that several important problems of the IHM in aerospace 
applications render themselves into a set of coupled stochastic differential 
equations ideally suited for analysis in the framework of dynamical inference. 
Specifically, we show that an open problem of the SHM for an on-board de- 
tection of the stage separation failure and damage in composite materials can 
be reduced to the DI of a set of oscillators representing the dynamics of the 
amplitudes of the eigenmodes. In fact, fast online DI of a set of oscillators 
with time varying parameters is of particular importance across a wide range 
of interdisciplinary applications including e.g. neurophysiology [30] and active 
control of combustion instability in liquid motors [7]. 

1.8.1 General inferential framework for a set of coupled os- 
cillators 

In this section, we briefly outline a general framework for a dynamical 
inference of such sets of oscillators. We consider a system of equations in the 
form system in the form 

= oijVj + bjVj + Cj v | - <ij + r/j + yl),, £j, 

Qj = -PQj +lj vj, (1-71) 

<&(*) &(*')> = Sij 6{t - t'), j = 1: L. 

The system (1.71) represents a dynamics of L nonlinear oscillators with co- 
ordinates q j and velocities Vj driven by random force with components 
which are mixed by a diffusion matrix with elements Djj . This form of equa- 
tions that can model a wide range of nonlinear phenomena (see e.g. [5]) and 
can be reduced, in particular, to a set of linear oscillators obtained in previous 
two sections un the context of online SHM problem. 

In practice, the measurements of the dynamical variables Vj(t) is a sepa- 
rate important problem in each specific application, in this example we use a 
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measurement model in a form of measurement matrix X , which is frequently 
found in theory of dynamical inference, 


Hi. — X \j Vj . 


(1.72) 


Here pi are measured variables, which are related to Vj by linear transfor- 
mation with an unknown matrix X. All the accessible information is con- 
tained in yi. The problem is, therefore, to learn model parameters A4 = 
{??,:, a,, gi(0), 7i, Dij,Xij} from time series data {yj. 

To develop a general inferential framework for a solution of this problem 
we notice that the second linear ordinary differential equation in Eqs. (1.71) 
connecting coordinates and velocities can be integrated explicitly (see [30, 16]) 
to obtain ^ 

qj{t) = 7 [ dTe~ l3 ^ t ~ T \j(T) + e~ l3t q j ( 0 ). 

Jo 

On substituting (1.73) into the first equation in (1.71) we have 

'■./ = a j v i + b j v j + c J v j + Vj 

- t j j dre-^^Vjfr) - 0 ) + 

Jo 

Here j = 1, ..., L and <?j(0) is a set of initial coordinates for unobservable vari- 
able qj ( t ) . Thus the reconstruction of unobservable variables q :t (t) is reduced 
to inference of L initial conditions qj{ 0). 

Variables Vj(t) can also be excluded from further consideration by using 
eq. (1-72). Indeed, on substituting v = X~ 1 y into (1-74) we obtain in vector 
notations: 


(1.73) 


(1.74) 


y = Xa (X- 1 y) + Xb (X~ x yf + Xc (X" 1 yf + e~ pt Xq Q 

-[ e /3(t_r) X 7 (X -1 y) dr + Xr) + xVDf(t), (1.75) 
Jo 


where a, b , and c are vectors with components {07}, { } , and {cy } respec- 
tively, qo = q(t = 0) and 
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— 1 % liVi 


0 \ 
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The Eq. (1.75) can be rewritten in a more general form 


Vi — &ij Vj T bijk Vj Vk T Cijkl Vj VkUl T ^]i T 
-$iiuyi - e~ /3t ; 


(1.76) 
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where the meaning of the parameters with ~ is clear from the comparison of 
(1.75) and (1.77), see also [30, 16]. 

Note that presentation (1.75) covers the whole model space of a set of 
coupled nonlinear oscillators with polynomial base functions of the power 3. 
This presentation can also be viewed as an independent model for inferring 
unknown parameters of a system of coupled nonlinear oscillators in the case 
when oscillators velocities yi can be measured directly. 

To infer the set of unknown parameters 

-'Xl — { !// , eXij , bij , C,j /,:/ , O'/j , Zj , Dij } 

of the system (1.77) of L coupled oscillators (1.75) within our inferential frame- 
work 1.2 one can introduce the following base functions 

<t>{x) = {i, yi, ...,y L ,yl, ym, -,yiyL,yl, V2V3, -,y2yL, -,yl, ^ 

y\y 2 , y\yL,vl, y\y 1 , •••, vIvl, •••, vIvl-i , e~ pt } 

where 4>j = f* yi(r)e^ T ~ t \ 

The total number of unknown parameters is N tot = (2 + 2 L + L 2 )L + 
L(L + l) 2 /2 and is increasing as L 3 with the dimension of the system. These 
parameters can be inferred directly from (1.77) using Eqs. (1.8) - (1.11) is the 
measurement matrix is known. In the case of unknown measurement matrix 
one has to infer additional L 2 coefficient of unknown matrix Xij , which is not 
symmetric in general. To solve this problem we notice that in practice the 
number of the coefficients of the original system is always significantly smaller 
then the full set of equations N to t, because of the symmetry, which is always 
present in real systems. For example, for the sets of linear oscillators that 
appear in eigemode expansion of the SHM problem in the on-board analysis 
of stage separation failure and damage of the composite materials the 2 nd 
equation in (1.71) has the form qj = Vj, i.e. (3 = 0, 7 j = 1, and all bijk and 
c-ijki are also zero. Furthermore, the eigenfrequencies of expansion in eigen- 
modes of every specific structural component that requires on-board SHM 
must be learned beforehand from extensive preliminary ground and in-flight 
tests, which significantly simplifies the problem. 

1.8.2 Numerical example 

We now consider a numerical example with a 2D set of oscillators (1.74) 
where all the coefficients of (1.74), but f3 are assumed unknown together with 
unknown coefficients of the measurement matrix {a+ }. In 2D case the set A4 
of variables of the transformed dynamics (1-75) corresponds to the following 
set of the base functions 

<i>(x) = {! , 2 / 1 , V 2 , yi, v\ , 2 / 12/2 , 2/1 , 2 / 1 2/2 , 2/12/1 , y\ > $ 1 , $2 , e~ m }. 

In general it should be possible to reconstruct all unknown coefficients of the 
original system for any number of FHN oscillators as long as we can establish 
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the connection between the set Xi and the set of original unknown variables 
of (1.71) 

•Ad — { Cf j . Oil > bi , Ci , , '/i - (Ji : (0 ) , Di j , X{j } 


Here we introduce explicit relations for the case L = 2. 


m 
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= X~ 1 
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DX" 1 = X -1 !}. (1.78) 
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(1.79) 

(1.80) 


The unknown elements Xij of the inverse measurement matrix X^ 1 and pa- 
rameters with tilde are the model parameters of the transformed system (1.75) 
that can be inferred directly using time series data {?/,}. Relations (1.78)-(??) 
allows one to reconstruct 15 unknown parameters of the original system, in- 
cluding elements of the noise and measurement matrixes. Note, however, that 
coefficients (1 + cq) can also be assumed unknown in general and the following 
relations can be used to reconstruct them 


1 + Oi\ 
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*11 

2 a.’n^i 2 

*12 

= x- x 

bin 
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bl 22 
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Similarly, the relation between the coefficients for polynomials of power 3 are 
given by 


*11 

2 xhxi 2 

2 xiiXi 2 

*12" 
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2X21X22 

2 x 21 X % 2 
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* 22 j 
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C212 
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C222 


(1.82) 


We now analyze the convergence of the method in the case when all the 
parameters of the reduced model (1.75), including elements of the measure- 
ment matrix are unknown. The sampling rate was 35 kHz, we used 9 blocks of 
data with 5000 points in each block, and these block of data were generated 
at random 1000 time to analyze the statistics of the convergence. 

To reconstruct both the mixing matrix X and the parameters of the orig- 
inal system Xi from the inferred parameters Xi of the transformed system 
(1.75) we have to solve equations (1.78-1.82) with respect to elements of M. 
In a particular case of transformation given in a simple form of Eqs. (1.78- 
1.82) the solution of this problem can be found using standard nonlinear least 
square method [8], however, an additional optimization over the set of initial 
values may be required. We stress that the present technique is not restricted 
to the 2D case and can be employed in the general case of L oscillators. 

The results of the inference of the transformed and original parameters 
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TABLE 1.5: Values of some of the transformed and original 

coefficients inferred using 30000 points obtained from measurement matrix 
and real parameters reconstruction. The actual values (second column) 
are compared with the inferred values (third column), relative errors are 
given in the last column. 


Parameter 

Actual Value 

Estimate 

Error 

Vi 

0.9200 

0.924384 

0.022624 

r P 

0.3500 

0.351001 

0.009063 

&222 

1.7550 

1.758011 

0.037047 

&112 

-2.1086 

-2.114731 

0.068268 

A n 

1.7 

1.686459 

0.796526 

V-|2 

0.8 

0.794263 

0.717092 

A 21 

0.2 

0.196746 

1.626811 

A 22 

0.9 

0.898222 

0.197610 

m 

0.4 

0.406227 

1.556788 

m 

0.3 

0.302462 

0.820660 


-0.35 

-0.351992 

0.569082 

OL 2 

-0.2 

-0.200376 

0.188228 

h 

1.35 

1.357427 

0.550145 

^2 

1.2 

1.203863 

0.321885 

Cl 

-1.0 

-0.999520 

0.047957 

C2 

-1.0 

-0.999114 

0.088582 


of a set of oscillators are summarized in Tab. ??. It can be seen from the 
table that the DI method allows us to reconstruct both the elements of the 
measurement matrix and the parameters of the transformed (1.77) and original 
(1.74) systems. Our analysis shows that the relative error of inference better 
then 2% is achieved in less then 1 sec of measurements. 

We, therefore, conclude that the general inferential framework for a set of 
coupled oscillators provided by Eqs. (1.8) - (1.11), (1.77) - (1.77) is one of 
the primary candidates for on-board SHM in aerospace applications. 


1.9 Conclusion 

It was shown that several important applications of IHM in aerospace ren- 
der themselves into a problem of on-line inference a set of coupled stochastic 
differential equations. A convenient approach to a solution of this problem was 
introduced within the framework of dynamical inference. It was shown that in 
a wide range of practical cases this problem can be solved analytically using 
path-integral approach to the maximum likelihood estimation, the robustness 
of the solution in the presence of strong dynamical noise was demonstrated in 
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application to the inference of archetypal nonlinear stochastic Lorenz system. 
The ability of the method to infer simultaneously discrete and continuous 
parameters was further illustrated by the application of the DI to the fast 
on-line estimation of parameters and detection of leaks in stochastic hybrid 
three tank system, which standard engineering control benchmark, the later 
problem was considered in a specific context of ground support system for 
liquid fuel filling system. 

In the remaining part of the chapter a progress in development of some se- 
lected FD&P systems for in-flight structural health monitoring was discussed. 
In the first example a dynamical model of the SRM ballistics in nominal and 
off-nominal regimes was presented. It was shown that the DI approach can be 
used to detect abrupt fault-induced changes in the parameters of the SRM. It 
was shown further that the DI method can be used to reduce the probability 
of “misses” in detection of the SRM faults on-board by inferring coefficients of 
the diverging terms in polynomial fit to small deviations of pressure from the 
nominal regime and generating an earlier warning when the inference results 
begin to converge. 

In the next example a progress in development of in-flight diagnostic and 
prognostic system for stage separation failure was considered. A mathematical 
formulation of the problem of structural dynamics of a truncated cone under 
unsymmetrical unsteady load due to impact was given that provides estimates 
for the stresses required for nozzle extension buckling. The results of the high- 
fidelity simulations in ABAQUS confirmed the buckling scenario and allowed 
to establish correlation between the strength of the impact induced torque and 
the damage of the nozzle. The model of the thrust vector control was used 
to predict the tilt/rock angles and forces generated by the TVC in response 
to impact. A dynamical analysis of these forces allows one to build FD&P 
system for stage separation failure and to cast the problem in terms of the DI 
of a set of driven oscillators representing amplitudes of the nozzle extension 
eigenmodes. 

In the third example a problem of in-flight structural health monitoring 
of composite sandwich plate was analysed. It was shown that normal mode 
expansion method is capable of capturing the physics of wave propagation 
in plate structures and may describe damage state of the plate. As a result, 
dynamic inference method is very promising for SHM of the plates. It was 
demonstrated in this context that it is sufficient to keep in the expansion up 
to ten terms to grasp main characteristic dynamical features corresponding to 
the damage present in the plate. 

Finally, the a general inferential framework for fast on-line reconstruction 
of the parameters of a set of coupled oscillators was presented in the context 
of the SHM problem. It was shown that the method allows for an accurate 
reconstruction of the system parameters even when the velocities of teh os- 
cillators are mixed by the linear measurement matrix with unknown elements 
and coordinates of the oscillators are not accessible for measurements. We 
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concluded that the DI is one of the primary candidates for on-board SHM in 
aerospace applications. 

Ultimately IHM research will yield integrated, multi-disciplinary analy- 
sis and optimization capabilities that enable system-level designs providing 
graceful recovery from in-flight failures, computationally efficient tools for in- 
flight prognosis of aircraft health including integrated predictive and sensor 
capabilities, and preventative and adaptive systems for in-flight operability 
and informed logistics and maintenance. 
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